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Combining incoming and outgoing characteristic formulations can provide numerical relativists 
with a natural implementation of Einstein's equations that better exploits the causal properties of 
the spacetime and gives access to both null infinity and the interior region simultaneously (assuming 
the foliation is free of caustics and crossovers). We discuss how this combination can be performed 
and illustrate its behavior in the Einstein-Klein-Gordon field in ID. 

on : 

Q\ ■ I. INTRODUCTION 

G\ : 

In recent years the characteristic formulation of G.R. jt]]^] has proven to be a valuable tool for numerical relativists^]. 
' Its use has enabled researchers to: study critical phenomena 4.5J; obtain the first 'unlimited' evolutions of generic 
single black holes pj; study accretion of matter by black holes |.9 (via a formulation based on foliating the spacetime 
with incoming nulTnypersurfaces) and obtain practical access to future null infinity where the gravitational radiation 
and other physically relevant quantities are obtained unambiguously [ [10[|llf| (via an outgoing null hypersurfaces 
foliation). These investigations adopted either an outgoing or an incoming nulTnypersurfaces foliation of the spacetime. 
Regrettably, in choosing one or the other one has access to future null infinity or to the black hole interior (as the 
y—i , incoming null surfaces can penetrate it but not the outgoing!) but not both. This, could be regarded as a drawback for 
• the formulation as it is desirable in many cases to access both regions; for instance, to describe a black hole accreting 
(*C) ' mass from its surrounding and the gravitational waves produced in the system. 

^C) ' At present, obtaining a model capable of addressing both problems at the same time requires using Cauchy- 
characteristic matching (CcM) — p~3| or the conformal formulation of Einstein's equations Jl^-bSt. Unfortunately, 
y—i . despite much progress in both approaches neither of them is fully working in 3D. The main problem facing the CcM 
■ approach is the involved set of interpolations back and forth from a 3D cartesian grid to 2-spheres at the matching 
j worldtube. On the other hand, the conformal approach has as a drawback its large number of variables which imposes 
, further constraints on present computer capabilities. 

The characteristic formulation |l,|], on the other hand, employs a minimum set of variables (6 compared with 12 in 
^ ■ the ADM equations or 53 in the conformal equations). The success obtained with it in both formulations suggests that 
a combination of both could well be a natural way to study both regions. This approach, which we call "characteristic- 
. characteristic matching" (c 2 M) divides the spacetime in two regions. The inner region (inside some worldtube T) to 
be studied by an incoming characteristic formulation and the outer region where and outgoing formulation will be 
used. Interpolations are carried out at T, however in 3D it only involves two-dimensional interpolations. In a related 
. £h \ approach, the double null formulation |17[] , Einstein equations are also integrated along incoming and outgoing null 
, hypersurfaces and its implementation might display similar stability properties as well. Unfortunately, to date, there 
are no numerical implementations of this approach in 3D but studies on how to carry out such a task are under 
development [Q. Naturally, neither of these approaches is expected to be applicable to all situations since, for 
instance, the presence of caustics prevents one to use them. However, in a caustic-free foliatiorjT], c 2 M represents 
an enticing proposal in light of the successful application of the characteristic formulation both in the vacuum and 
matter cases. 

We briefly describe the characteristic formulation (both in the incoming and outgoing formulation) in section 2. 
Section 3 outlines the matching strategy to be used in this work. We describe the numerical implementation in section 
4 and test it in section 5. Finally we conclude with some comments and future research directions. 

II. THE CHARACTERISTIC FORMULATION OF G.R. 

Einstein's equations can be expressed in notational form as G a b = S>irT a b. Where G a b is the Einstein tensor and T a b 
the stress e nerg y tensor of the matter distribution. In the particular case of a massless scalar field $ coupled to G.R., 
T ab results ffj 



*For a recent review on its use in numerical relativity see j|. 

t Aside from the natural caustic present at the origin of the null cones used to coordinatize the incoming null surfaces. 
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T ab = V a $V h $ - l/2 5ab V c $V c $ . 



(2.1) 



We introduce a coordinate system adapted to either incoming or outgoing null hypersurfaces in the following 
way: the outgoing (incoming) lightlike hypersurfaces are labeled with a parameter u (u); each null ray on a specific 
hypersurface is labeled with x A (A — 2, 3) and r is introduced as a surface area coordinate (i.e. surfaces at r — const 



have area 47rr 2 ) (see Fig. 1). In the resulting x a 
§1 



(u(v),', 



coordinates, the metric takes the Bondi-Sachs form 



ds 1 = - 



3 V/r - r 2 h AB U A U B ) du 2 - 2e 2f3 dudr 



-2r 2 h AB U B dudx A + r 2 h AB dx A dx B , 
for the outgoing case, while for the incoming case: 



ds z 



where in both line elements, h AB h B c 
incoming line element, as discussed in J 
This fact was used in Jgjy. to obtain a > 
for the outgoing case |M. 



r 2 h AB U A U B ) dv 2 



2e 2/3 dudr 



(e 2l3 V/r - 

-2r 2 h AB U B dudx A + r 2 h AB dx A dx 1 



(2.2) 



(2.3) 



S A and det(h AB ) — det(q AB ), with q AB a unit sphere metric. Note that the 
J, can be obtained from the outgoing version by the substitution (3 — > f3+iir/2. 
J implementation of the incoming null formulation from the one constructed 



null 
infiaiix 




/ N 



FIG. 1. Incoming (A) and outgoing (B) formulations. In (A), the interior of F is covered by a sequence of incoming light 
cones. In (B), the exterior of F is covered by a sequence of outgoing light cones. 

In the spherically symmetric case, Einstein's equations and the evolution equation for the scalar field reduce to: 



/3 r = 2vrr$ 2 r . 



Vr 



2(r$), 

for the outgoing case, while for the incoming one, 



r 



Vr 



=,2/3 



2(r$) lPO = -(r7$, r ), r 
r 



(2.4) 
(2.5) 

(2.6) 

(2.7) 
(2.8) 
(2.9) 



(where we distinguish with a tilde the incoming fields from the outgoing ones). 

Given that initial data for <I> ($) on an initial outgoing (incoming) hypersurface and at some worldtube T consistent 
values of V, (3 and $ are provided, the equations can be straightforwardly integrated to yield a unique solution 
p9|,plj. Consistent boundary data on T must satisfy constraint equations obtained from G r a = 8ttT^ (a — u,x A ). 
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These conditions are automatically satisfied when matching across T to an interior (exterior) solution generated by 
the incoming (outgoing) evolution |10| . 

Asymptotic quantities of special interest are the Bondi mass and the scalar news function p2| 

M(u) = ^e- 2H r 2 {V/r). r \ r=x , (2.10) 
N(u) = e- 2H r$, u | r=00 ; (2.11) 
where H = f3\ r=00 . The Bondi mass loss equation is 



e 



- 2H M, U = -AttN 2 (2.12) 



This last relation provides an important test to the accuracy of our evolution as 



M(u = u ) = M(u = u f ) + 4vr / e 2H N 2 du , (2.13) 

J u 

which we monitor to assess the validity of our numerical implementation. 

III. MATCHING 

Our aim is to match both formulations at a worldtube T defined by r = R = const (for simplicity, in fact this 
condition can be relaxed to r = R(u) and in the 3D case to R{u,x A )). A n incoming formulation to be used in the 
region d — {(v,r) 7 r G [0, Ri]} and an outgoing one in the region C = {(u,r),r e [i?2,oo)} (R2 < R < Ri). By 
coupling these codes, consistent boundary values for the outgoing case can be obtained from the incoming variable 
via a coordinate transformation and viceversa. In the spherical case, obtaining this coordinate transformation is 
straightforward^]. Recall that the incoming line element is given by 

d s 2 = e 2 ^dv 2 + 2e 2 ^dvdr + f 2 dfl (3.1) 



Clearly, the transformation 



dr = df (3.2) 
du = dv — V jfdf (3-3) 

takes the incoming formulation into the outgoing one. This transformation is valid everywhere and can be used to 
obtain boundary values by simply using the transformation law for a 2-tensor field obtaining: 

= e /3 , V = -V . (3.4) 

Additionally, the scalar field must be continuous across T, so the condition $r = applies. Also, in order to ensure 
"sheet sources" produced by discontinuities in derivatives across T will be not be present one might further impose 

(fc Q V Q $)- = (fc a V a $)+; (3.5) 

where the superscripts — and + denote the regions to the left and right of T respectively and k a is an arbitrary vector 
at the worldtube boundary with non zero radial component. 

A. Singularity Excision 

In order to avoid dealing with singularities, we use excision to remove it from our computational domain. This 
technique first suggested by Unruh fl23]l is at present the most successful approach to handle singularities. As it is 
customary, we excise from the computational domain the region inside an inner trapping boundary defined by the 



The matching strategy in the 3D case is naturally more involved and we comment how it can be implemented in the appendix. 
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following. Given a slice S on an incoming null hypersurface Af described in null coordinates by r = R(v,x A ), the 
divergence of the outgoing null normals is given by 

= —V - -L{^h AB R B - r 2 U A )], A 

-r{r~ l e 2p h AB ), r R A R B + r 2 U A R A . (3.6) 

Hence, the slice is marginally trapped if 6/ = 0. In the spherically symmetric case, this reduces to r = R(v), and the 
divergence of the outgoing null rays results 

e -2/3 

F^-2V—, (3.7) 

Therefore, in our implementation we simply look for the largest r — r e for which F < 0; for r < r e the field variables 
are not integrated and simply set to the values they have at r e . This straightforward procedure excises the singularity 
from the computational domain. 



IV. NUMERICAL IMPLEMENTATION 



We introduce a compactificd radial coordinate x = r/(l + r) (so that x = 1 — > r = oo). Defining the x location of 
the worldtube by x wt = -R/(l + R), we introduce the outgoing radial grid as X0i = x wt + (i — 2)Axo where Axo — 
(1 — x wt )/ (N xo — 2). Analogously, we define the incoming radial grid as xij — (j — l)Ax, where Axi — x wt / (N x i ~ 2). 
Notice that XO2 = xi^ !ci _ 1 — x wt - Hence, there is an overlapping region which will be used to interpolate from one 
formulation to the other in order to define boundary values. Additionally, we choose v — u at T. 



A. Evolution Scheme 



Integration of the evolution equations is done via the parallelogram approach introduced in p4|. We here describe 
briefly what this method entails for the outgoing case (it is str aigh tforward to translate it to the incoming case so we 
do not expand on it). Basically, one recognizes that equation ( |2.6| ) is equivalent to 

2(r$), ur - (-(r*), r ) =-(-) *; (4.1) 



t J „ \ r 



defining G = r$, the previous equation can be rewritten as 



□G=-(T) -■ (4-2) 



r J „ r 



This equation corresponds to the 2-dimensional wave equation in the (u, r) plane and can be easily integrated over 
a null parallelogram defined by two sets of outgoing and incoming null geodesies (see figure ^ . Hence, the integral 
identity 

G Q =G P + G R -G r + \ f dudr(-), r - , (4.3) 

2 J a r r 

can be used to integrate the scalar field up to global second order accuracy as described in p^j . Additionally, the 
"hypersurface" equations are discretized by centered second order differences yielding a global second order accurate 
numerical implementation. 
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FIG. 2. Parallelogram algorithm set up. The null parallelogram is defined by the intersection of 2 outgoing geodesies with 
2 incoming ones. 



B. Matching scheme 

A straightforward matching implementation can be obtained by the following scheme. Suppose that at time v = 
u = u a (= nAu) all field variables are known (refer to figure 3). In order to obtain starting values at the level 
v = u = u + Au one can find the intersection of the incoming null surface at v = u + Au with the outgoing null 
surface corresponding to u = u , this intersection can be easily obtained in the following way. First, consider the 
outgoing null ray emanating from x Q = (u a , x wt ) (R in figure 3); the intersection point Xi of this ray with the incoming 
one going through at u — u a (indicated by S in figure 3) is obtained by requiring that 



ds z = 



(4.4) 



which implies that dr Q — ^du. Hence, the values of the fields at Xi can be easily obtained by interpolation of the field 
on the outgoing patch. This value will be used as starting values to carry out the incoming radial integration. By 
proceeding in an analogous way one can get values at x" — (u + Au, R — dr ) (P in figure ffl) which provide starting 
values for the radial integration on the outgoing null surface. Additional ly, starting values for the hypersurfaces 
equations are also obtained by interpolations and applying transformation (3.4). 

To ensure continuity of the field and its radial derivatives we proceed as follows. Since the intersection of the null 
surfaces naturally forms a null parallelogram (see figure 3), one can integrate the evolution equation for $ (as in 

section 4.1) and obtain the value of $2 +1 (since the values of $ at all the other corners are either known directly or 
can be obtained by interpolation from the values of the fields at the previous level). This procedure naturally ensures 
continuity of the field and its derivative across T. 




FIG. 3. Matching the characteristic formulations. Black squares indicate grid points while the grey circles indicate the 
intersection of incoming and outgoing null rays. Points P, Q, R and S define the corners of the parallelogram that will be used 
for the integration to obtain the value of $ at the n + 1 level on V. 
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V. TESTS 



Our purpose is to present a model demonstrating the feasibility of a stable algorithm based upon two regions that cover 
the spacetime, where equations are expressed with respect to and incoming or outgoing characteristic formulation. 
The inner region is evolved using an incoming null algorithm whose inner boundary lies at the apparent horizon (in 
case it exists, otherwise at r = 0) and whose outer boundary R\ lies beyond the inner boundary of the outer region 
(located at R2) handled by an outgoing null evolution. 

A. Scalar wave on a flat background 

As a first test, we choose initial data corresponding to a scalar wave on a flat background. The initial data for the 
scalar wave is described by 

fi = 0; V = -r; $ = 0; (5.1) 

$ f a(l - r/R a )%l - r/R b f for r G [R a , R b ] (5 2) 

\ otherwise ^ ' ' 

where the values of (3 and V are determined by matching at R = 5 and integrating the hypersurface equations. 
The value of the field at the origin is chosen so that <f> jr = 0. The parameter a was set to be small enough so that 
the wave would reach the origin and be dispersed away without a black hole being formedPL In our runs, we choose 
R a = 10, Rb = 20 and study the system for differ ent values of a (with N X i — N xo = 162). As a test of our numerical 
implementation we checked that equation ( 2.1 3|) was satisfied throughout the evolution. Our results for the case 



a = 1 , are illustrated in figures [| and ^. Energy conservation is satisfied and the evolution of the scalar field proceeds 
as expected. The pulse is originally registered at the outgoing patch; as it travels towards the origin it crosses T, 
reaches the origin and it disperses to infinity radiating all the energy contained in the initial data. The spacetime at 
late times is flat and all the energy of the initial pulse has been radiated away. 



p 

P+<W„ 



FIG. 4. Conservation of Energy. The plot shows the value of the Bondi mass, the total radiated energy and the addition of 
the Bondi mass at a given time with the energy radiated until that time (indicated by the solid line) . Clearly, the total energy 
is conserved throughout the evolution. 



§ In principle, we could use our numerical implementation to try to look for critical phenomena as first observed by Choptuik 
25l, we defer such a study to a future work. 
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FIG. 5. A sequence of snapshots of the field on both the interior (left) and exterior region (right). From top to bottom the 
snapshots correspond to u — 0, 17, 34 and 51. Initially the scalar field is non zero only on the exterior region, at later times the 
field crosses V, reaches the origin and it is radiated away to future null infinity leaving behind a flat spacetime. 



B. Collapse of a spherically symmetric scalar wave onto a black hole 



We choose initial data corresponding to a Schwarzschild black hole of mass M which is well separated from a localized 
pulse of (mostly) incoming scalar radiation. This problem has been thoroughly studied in the past by Marsa and 
Choptuik |2q| in the Cauchy formulation of G.R. Here, we do not try to extend their results, rather, to verify the 
applicability of c 2 M . Our configuration is such that, initially, there is no scalar field present on the incoming null 
patch, so the initial data are simply, 

l> = 0,/3 = (5.3) 
V = 2M- r. (5.4) 



On the outgoing patch we choose $ as a pulse with compact support given by (5.2). The values of (3 and V are 
determined by matching at R and integrating the hypersurface equations. Inner boundary conditions are not needed 
since the region inside the marginally trapped surface is excised from the computational domain. 

In our tests, we choose R = 20M, R a = 40M, R^ = 8QM and again study the evolution of the system for different 
values of a (with N X i = N xo = 162). As an illustration, we present the results obtained with a = 2. As shown in 
figures || and |?], the scalar field propagates towards the origin, part of its "original" energy falls into the hole which 
increases its mass and the rest is radiated away. 



0.6 - 

M B 

— p 

0.4 - P+M. 



0.0 20.0 40.0 60.0 80.0 100.0 



FIG. 6. Conservation of Energy II. The plot shows the value of the Bondi mass, the total radiated energy and the addition 
of the Bondi mass at a given time with the energy radiated until that time (indicated by the solid line). Values for this run 
were a = 2, R a = 40M, Rt = 80M and R = 20. Part of the initial energy contained in the pulse is radiated away while part 
of it is accreted by the black hole. The final spacetime corresponds to a Schwarzschild spacetime with total mass larger than 
the initial hole. 
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1.2 - 




0.4 I ' ' ' ' ' ' ' ' ' 1 

0.0 20.0 40.0 60.0 80.0 100.0 

U 

FIG. 7. Apparent horizon location vs. time. The figure shows how the apparent horizon location grows as energy is 
accreted by the hole. The final asymptotic value of 1.1381 agrees well with twice the value of the Bondi mass at late times 
(2Mb = 1.1388). 

The code successfully excises the hole and the area of the apparent horizon reaches an asymptotic value which 
agrees with twice the value of the Bondi mass at late times (see figure 0). This result fully agrees with what is 
expected, since at late times the spacetime must be describable by aSchwarzschild metric as a consequence of the 
no-hair theorems. A sequence of this evolution is illustrated in figure ||. 



0.14 — 
0.06 - 
-0.02 — 




X X 

FIG. 8. A sequence of snapshots of the field on both the interior (left) and exterior region (right). From top to bottom the 
snapshots correspond to u — 0, 15, 30 and 45. Initially the scalar field is non zero only on the exterior region, at later times 
the field crosses F, and reached the black hole which accretes some of the energy and so the masked region increases (shown 
by the constant lines in the second to fourth plots on the left). The rest is radiated away to future null infinity leaving behind 
a Schwarzschild spacetime with a larger mass. 

Another important test of the algorithm focuses on general aspects of wave propagation in a curved background. 
The qualitative features of the evolution are well known p7| , p8[ . In particular, a distinctive aspect of wave propagation 
on curve backgrounds is that Huvgens principle does not apply. In a non-flat spacetime, radiation backscatters giving 
rise to quasinormal modes and tails [£9|. Notably, for a given Schwarzschild potential, the power of the late-time tail 
is characteristic of generic compact support initial data and depends on the multipole structure of it. In particular, 
at future null infinity the dependence on time goes as u~ 2 ~ l (with 1 the multipole index). In our case, I = since 
the initial data corresponds to a spherically symmetric pulse; the value measured from the simulation is —2.1 which 
agrees with the expected value of —2 (see figur e B [) . The results obtained here reproduce and, in a way, expand those 
obtained with Cauchy-characteristic matching |30| , as in this case we dynamically find the horizon and use it to excise 
the hole. The behavior of the apparent horizonlocation gives a good indication of the energy absorbed by the black 
hole and its area at late times agrees with that obtained by the asymptotic value of the Bondi mass. 
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FIG. 9. Late-time power-law tail at null infinity. The measured slope of —2.1 compares well with the theoretical value of 



VI. PRESENT AND FUTURE OUTLOOK 



This work shows that the matching approach of two characteristic formulations provides a way of treating both the 
inner region and the outer region. At the inner region, it tackles the excision problem as well as other treatments 
like Cauchy characteristic matching and purely Cauchy or characteristic treatments. It has some advantages over the 
Cauchy approach because its implementation is much simpler and computationally cheaper. Its disadvantage though, 
is that it can not be applied in the case where caustics might be present in the integration region. Additionally, 
it provides a means to reach future null infinity, something that only can be done with the help of a characteristic 
formulation or the conformal formulation of Einstein's equations. 

This study is the first in a line of research in numerical relativity aiming to fully match characteristic codes to study 
matter coupled to G.R. Preliminary works have shown that, where applicable, the use of a characteristic evolution 
can greatly help in investigating this problem. To date however, these studies have been restricted to either the 
incoming or outgoing formulation, thus not being able to study the inner region or obtaining physical quantities at 
future null infinity simultaneously. A combination of both approaches appears to be an ideal candidate since it would 
combine the main advantages of both directions giving access both to future null infinity and to the horizon region 
in a natural way. Additionally, it might help "avoid" caustics that can be encountered using the single mapping of 
each formulation. Just as several coordinates patches are necessary to describe nontrivial topologies, a combination 
of outgoing and incoming coordinates can help to obtain a smooth mapping of the manifold under study. The 
present work represents the first step towards achieving such a combination. At present, 3D robust implementations 
exist for both the incoming and outgoing characteristic formulations 10 fj]; further, matter treatments studied with 
these formulations have been successfully implemented h|. Thus, demonstrating the feasibility of combining both 
approaches in ID inspired this work. Future studies wilTaddress the matter case in ID (ie. implementing the fluid 
equations rather than that of a massless scalar field) and proceed with full matching in 3D. 
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VIII. APPENDIX 



A matching strategy for the general 31? case can be devised in a way similar to the one implemented in the Cauchy- 
characteristic matching approach. In the c 2 M case, however, its application is simpler since no interpolation from 
a 3D cartesian grid to 2D spheres (slices of the matching worldtube) is required. This matching technique [|l3| is 
straightforward. First assume that all field variables are known at hypersurfaces Af v and Af u -Au (and earlier ones); 
that the matching worldtube is located at f = R = const and that u = v at T. The first step which is referred 
to as extraction involves obtaining start-up values (field values at the worldtube) for the outgoing variables at Af u . 
These values are obtained by transforming the incoming metric tensor at Af v on T to the outgoing coordinate system. 
Instead of solving the geodesic equation to obtain this coordinate system one can proceed as follows. Choose an 
"affine" coordinate system x a — (u,X,x a ) by 







x a \ r + Xl a + 0(A 2 ) 



(8.1) 



where x a |r 
is at r = R 



(v,r = R,x a ); l a is an outgoing null vector (normalized by l a dv a — —1) and we assume the worldtube 
const. By inspection, equation 



.1 



) is the solution to the geodesic equation with affine parameter A 
up to second order in A. Using this transformation, one can obtain the metric in outgoing null affine coordinates by 
a simple tensor transformation. The final step requires obtainin g r. as its value can only be known after the angular 
components of the outgoing null metric are obtained, (refer to ]13| for further details). With these values, one can 
integrate Einstein's equation and completely determine the field values on 7V M . 

The second step, referred to as injection requires us to obtain the field values at the i nter section of A/^+ad with 



J\f u . This intersection can be obtained up to second order accuracy by means of equation (S.l), since 



v + Av = v + \l°\ r + 0{\ 2 ) , 
n = R + Xl r \ r + 0(X 2 ). 



(8.2) 
(8.3) 



From which we can deduce the value of r%. Thus, since the values of the fields are known everywhere on M u , we can 
i nter polate to obtain the values of the outgoing fields at the intersection point. Finally, by means of the transformation 
flS.ip , the values of the incoming fields are obtained with which one can integrate the equations inward. 
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